from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
from ngsolve.bem import *

keys: Acoustic Scattering, Brakhage-Werner formulation, Helmholtz equation, Dirichlet data

Helmholtz Problem

We consider a Helmholtz problem with a plane wave as initial condition, i.e., \(m = e^{i\,\kappa \,x} \) and

\[\begin{split} \left\{ \begin{array}{rcl l} \Delta u + \kappa^2 \, u &=&0, \quad &\Omega \subset \mathbb R^3\,,\\ \gamma_0\, u &=& m, \quad &\Gamma = \partial \Omega\, .\end{array} \right. \end{split}\]

Combined field integral equations combine single and double layer integral operators, one simple option is the Brakhage-Werner formulation.

The solution is represented as

\[ u = i \,\kappa \,\mathrm{SL}(j) - \mathrm{DL}(j)\,, \]

where \(j\) solve the boundary integral equation $\( \left( \frac{1}{2} + \mathrm K + i \, \kappa \,\mathrm V \right) \, \mathrm j = \mathrm u_{in} \qquad \text{on} \, \Gamma = \partial \Omega \)$

The single layer and double layer operator are defined using the fundamental solution of the Helmholtz operator,

\[ G(x,y) = \frac{e^{i\kappa \|x-y\|}}{4\pi\,\|x-y\|}. \]

The layer potentials \( \mathrm SL\) and \(\mathrm DL\) have the same form as in the Laplace case, except that the Laplace kernel is replaced by the Helmholtz kernel \(G\).

kappa=10
order=4

screen = WorkPlane(Axes( (0,0,0), Z, X)).RectangleC(15,15).Face()

sp = Fuse(Sphere( (0,0,0), pi).faces)
screen.faces.name="screen"
sp.faces.name="sphere"
shape = Compound([screen,sp])

mesh = shape.GenerateMesh(maxh=5/kappa).Curve(order)
Draw (mesh);

fes_sphere = Compress(SurfaceL2(mesh, 
                                order=order, 
                                complex=True, 
                                definedon=mesh.Boundaries("sphere")))
u,v = fes_sphere.TnT()
fes_screen = Compress(SurfaceL2(mesh, order=order, 
                                dual_mapping=True, 
                                complex=True, 
                                definedon=mesh.Boundaries("screen")))
print ("ndof_sphere = ", fes_sphere.ndof, "ndof_screen =", fes_screen.ndof)

ndof_sphere =  17010 ndof_screen = 31110

with TaskManager():
    C = HelmholtzCF(u*ds("sphere"), kappa)*v*ds
    u,v  = fes_sphere.TnT()
    Id = BilinearForm(u*v*ds).Assemble()

lhs = 0.5 * Id.mat + C.mat
m = exp(1j * kappa * x) 
rhs = LinearForm(m*v*ds).Assemble()

gfu = GridFunction(fes_sphere)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
with TaskManager():
    gfu.vec[:] = solvers.GMRes(A=lhs, b=rhs.vec, pre=pre, maxsteps=40, tol=1e-8)

GMRes iteration 1, residual = 60.52440579389444     

GMRes iteration 2, residual = 24.00171041505929     

GMRes iteration 3, residual = 10.816186639702561     

GMRes iteration 4, residual = 5.178501301182413     

GMRes iteration 5, residual = 2.549268448878952     

GMRes iteration 6, residual = 1.2888054956882669     

GMRes iteration 7, residual = 0.6549832602491108     

GMRes iteration 8, residual = 0.3554164558330167     

GMRes iteration 9, residual = 0.19990866606606375     

GMRes iteration 10, residual = 0.10836869993844714     

GMRes iteration 11, residual = 0.05931445314039659     

GMRes iteration 12, residual = 0.03294608123183362     

GMRes iteration 13, residual = 0.01381913910824682     

GMRes iteration 14, residual = 0.007360642470274796     

GMRes iteration 15, residual = 0.0033138105281745656     

GMRes iteration 16, residual = 0.0019472886266694597     

GMRes iteration 17, residual = 0.0011529608192671487     

GMRes iteration 18, residual = 0.0006865451125338563     

GMRes iteration 19, residual = 0.00040647668374399725     

GMRes iteration 20, residual = 0.00023525302450811214     

GMRes iteration 21, residual = 0.0001256833512855774     

GMRes iteration 22, residual = 7.272844838705902e-05     

GMRes iteration 23, residual = 4.282868528421695e-05     

GMRes iteration 24, residual = 2.238947095910066e-05     

GMRes iteration 25, residual = 1.132365408107729e-05     

GMRes iteration 26, residual = 5.681772319703246e-06     

GMRes iteration 27, residual = 2.8154812159697885e-06     

GMRes iteration 28, residual = 1.3918984915180512e-06     

GMRes iteration 29, residual = 7.931632776571544e-07     

GMRes iteration 30, residual = 4.523061248529777e-07     

GMRes iteration 31, residual = 2.723734968721736e-07     

GMRes iteration 32, residual = 1.609611766851694e-07     

GMRes iteration 33, residual = 9.476709863742893e-08     

GMRes iteration 34, residual = 5.2145026155061434e-08     

GMRes iteration 35, residual = 2.970933842689033e-08     

GMRes iteration 36, residual = 1.7339940734545346e-08     

GMRes iteration 37, residual = 9.167410459416053e-09     

Draw (gfu, order=5, min=-1, max=1);

Postprocessing on screen

uscat = GridFunction(fes_screen)
with TaskManager():
    uscat.Set(1j*kappa*HelmholtzSL(u*ds("sphere"),kappa)(gfu) - HelmholtzDL(u*ds("sphere"), kappa)(gfu), \
              definedon=mesh.Boundaries("screen"))

print ("Scattered field")
Draw (uscat, mesh, min=-1,max=1, animate_complex=True, order=4);

Scattered field

uin = mesh.BoundaryCF( {"screen": m }, default=0)
print ("Total field")
Draw (uin+uscat, mesh, min=-1,max=1, animate_complex=True, order=4);

Total field